{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 初始化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import matplotlib.pyplot as plt\n",
    "# macOS\n",
    "# plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']\n",
    "# Windows\n",
    "plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']\n",
    "# 正常显示负号\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import rc\n",
    "from matplotlib import dates\n",
    "rc('mathtext', default='regular')\n",
    "# 显示所有变量/最后一个变量\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = 'all'\n",
    "# InteractiveShell.ast_node_interactivity = 'last'\n",
    "import pandas as pd\n",
    "pd.set_option('display.max_rows', None)\n",
    "pd.set_option('display.max_columns', None)\n",
    "import numpy as np\n",
    "# conda install scikit-learn\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from keras.models import Sequential\n",
    "from keras.layers import LSTM, Dense, Activation,Dropout\n",
    "from scipy import stats\n",
    "from pylab import *\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import silhouette_score, silhouette_samples, r2_score, accuracy_score, explained_variance_score\n",
    "from sklearn.cluster import KMeans\n",
    "from keras.models import load_model\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from xgboost import XGBRegressor\n",
    "from sklearn.linear_model import LinearRegression\n",
    "# from fitter import Fitter\n",
    "from sklearn import model_selection, metrics\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "import xgboost as xgb\n",
    "from sklearn.model_selection import train_test_split\n",
    "from xgboost import XGBRegressor\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import numpy as np\n",
    "from sklearn import svm\n",
    "import datetime\n",
    "from time import time\n",
    "# pip install EMD-signal\n",
    "# from PyEMD import EEMD\n",
    "# from minepy import MINE\n",
    "from sklearn.neural_network import MLPRegressor\n",
    "\n",
    "args_wea = ['year', 'month', 'day', 'hour', 'minutes', 'temperature', 'cloud', 'sun', 'uv', 'water', 'wet']\n",
    "\n",
    "pv_time_index = ['2018/1/2', '2018/1/22', '2018/2/10', '2018/3/1', '2018/3/20', '2018/4/10']\n",
    "load_time_index = ['2018/1/2', '2018/1/22', '2018/2/10', '2018/3/1', '2018/3/20', '2018/4/10',\n",
    "    '2018/5/15', '2018/6/20', '2018/10/1']\n",
    "time_str1 = '20220101'\n",
    "time_str2 = '20220430'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据读取"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 负荷\n",
    "load_data = pd.read_csv('./data/load_2018.csv', encoding='gbk')\n",
    "load_data = load_data.set_index('时间')\n",
    "load_data.head()\n",
    "data = np.array(load_data['实际功率'])\n",
    "data = data.reshape(len(data) // 96, 96)\n",
    "data = pd.DataFrame(data).unstack().reset_index()\n",
    "data.columns = ['X', 'Y', 'Z']\n",
    "data['X'] = pd.Categorical(data['X'])\n",
    "data['X'] = data['X'].cat.codes\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "surf = ax.plot_trisurf(data['X'], data['Y'], data['Z'], cmap=plt.cm.viridis, linewidth=0.2)\n",
    "_ = ax.set_xlabel('时刻', fontsize = 15)\n",
    "_ = ax.set_ylabel('日期', fontsize = 15)\n",
    "_ = ax.set_zlabel('负荷', fontsize = 15)\n",
    "_ = fig.colorbar(surf, shrink=0.5, aspect=10)\n",
    "# _ = plt.savefig(work_dir + '/../Image/负荷曲线3D.svg', format='svg')\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 负荷预测\n",
    "## 特征选取"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = load_data\n",
    "_ = plt.figure(figsize=(16,9))\n",
    "mms = MinMaxScaler(feature_range=(0, 1))\n",
    "# _ = plt.figure(figsize=(16, 12))\n",
    "mine = MINE(alpha=0.6, c=15, est=\"mic_approx\")\n",
    "count = 0\n",
    "for i, time in enumerate(load_time_index[:8]):\n",
    "    # 归一化\n",
    "    # 定性分析\n",
    "    a = mms.fit_transform(data[data.index==time]['实际功率'].values.reshape(-1,1))\n",
    "    b = mms.fit_transform(data[data.index==time][data.columns[i]].values.reshape(-1,1))\n",
    "    if i not in [0, 1, 7]:\n",
    "        count += 1\n",
    "        _ = plt.subplot(3, 2, count)\n",
    "        _ = plt.plot(a, color='b')\n",
    "        _ = plt.plot(b, color='r')\n",
    "        _ = plt.xlabel('时刻')\n",
    "        _ = plt.ylabel('功率/p.u.')\n",
    "        _ = plt.title(f'负荷与{data.columns[i]}')\n",
    "        # 定量分析，计算最大互信息数\n",
    "        mine.compute_score(a.flatten(), b.flatten())\n",
    "        print(f'实际功率与{str(data.columns[i])}MIC：{str(mine.mic())}')\n",
    "_ = plt.tight_layout()\n",
    "_ = plt.show()\n",
    "\n",
    "columns = data.columns\n",
    "index = data.index\n",
    "data = mms.fit_transform(data)\n",
    "data = pd.DataFrame(data, columns=columns, index=index)\n",
    "data.head()\n",
    "# 计算任意一日'p'与前2天至4天前3个小时12个时刻的MIC\n",
    "# 前4个时刻MIC大于0.8\n",
    "for i in range(2, 5):\n",
    "    m1 = data[96*20:96*21]['实际功率']\n",
    "    for j in range(4):\n",
    "        m2 = data[(96*(20-i)-j):(96*(21-i)-j)]['实际功率']\n",
    "        mine.compute_score(m1, m2)\n",
    "        print(f'负荷与前{i}天{j}时刻前负荷的MIC：{mine.mic()}')\n",
    "\n",
    "# 计算一周各因素的MIC\n",
    "for col in data.columns:\n",
    "    mine.compute_score(data['实际功率'][96*31:96*38], data[col][96*31:96*38])\n",
    "    print(f'负荷与{col}的MIC：{mine.mic()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 训练数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gra import gra_df\n",
    "\n",
    "def get_train_data_load(df, x_par_names, y_par_name):\n",
    "    '''\n",
    "    df: Dataframe类型\n",
    "    df_pre: Dataframe类型, 待预测日数据\n",
    "    x_par_names: list类型, 表中用于作为x输入的列名\n",
    "    y_par_name: list类型, 表中用于作为y输出的列名\n",
    "    '''\n",
    "    # 待预测日即为data中最后一天\n",
    "    x_valid = np.array(df[x_par_names][-96:]).flatten()\n",
    "    y_valid = np.array(df[y_par_name][-96:]).flatten()\n",
    "    gra_0 = {'x':x_valid, 'y':y_valid}\n",
    "    days = int(len(df) // 96)\n",
    "    # 根据灰色关联系数生成相似日样本集\n",
    "    lis = []\n",
    "    for i in range(days-2):\n",
    "        x = np.array(df[x_par_names][-(i+2)*96:-(i+1)*96]).flatten()\n",
    "        y = np.array(df[y_par_name][-(i+2)*96:-(i+1)*96]).flatten()\n",
    "        gra_1 = {'x':x, 'y':y, 'day':i+1}\n",
    "        # 两列数求灰色关联系数\n",
    "        da = pd.DataFrame({'origin':gra_0['x'], 'new':gra_1['x']})\n",
    "        gra = gra_df(da).iloc[0,1]\n",
    "        dic = {'gra':gra, 'x':gra_1['x'], 'y':gra_1['y'], 'day':gra_1['day']}\n",
    "        lis.append(dic)\n",
    "    # 根据gra值进行排序，取相似度最高的n天作为训练数据\n",
    "    data = sorted(lis, key=lambda x:x['gra'])[-5:]\n",
    "    similar_days = []\n",
    "    [similar_days.append(d['day']) for d in data]\n",
    "    print(f'相似日:{similar_days}')\n",
    "\n",
    "    x_train, y_train = [], []\n",
    "    for i in range(len(data)):\n",
    "        x_train.extend(data[i]['x'][j] for j in range(len(data[i]['x'])))\n",
    "        y_train.extend(data[i]['y'][j] for j in range(len(data[i]['y'])))\n",
    "    args_num = len(x_par_names)\n",
    "\n",
    "    x_train = np.array(x_train).reshape(len(x_train) // args_num, args_num)\n",
    "    x_valid = np.array(x_valid).reshape(len(x_valid) // args_num, args_num)\n",
    "    print(f'x_train shape{shape(x_train)}')\n",
    "    print(f'x_valid shape{shape(x_valid)}')\n",
    "\n",
    "    # 相似日负荷曲线绘制\n",
    "    _ = mms.fit_transform(load_data[['实际功率']])\n",
    "    for i in range(x_train.shape[0]//96):\n",
    "        xs = [x_train[i*96+j][0] for j in range(96)]\n",
    "        xs = mms.inverse_transform(np.array(xs).reshape(-1,1))\n",
    "        _ = plt.plot(xs, 'b')\n",
    "    xs = [x_valid[i][0] for i in range(96)]\n",
    "    xs = mms.inverse_transform(np.array(xs).reshape(-1,1))\n",
    "    _ = plt.xlabel(\"time\")\n",
    "    _ = plt.ylabel(\"power/kW\")\n",
    "    _ = plt.plot(xs, 'r')\n",
    "    _ = plt.show()\n",
    "\n",
    "    return x_train, y_train, x_valid, y_valid\n",
    "\n",
    "mms = MinMaxScaler(feature_range=(0, 1))\n",
    "data = load_data\n",
    "# data = data.drop({'temperature', 'cloud', 'sun', 'uv', 'water', 'wet'}, axis='columns')\n",
    "data = data.drop({'时节', '月份', '风速', '降雨'}, axis='columns')\n",
    "# 由确定的最优特征集生成训练数据集\n",
    "x_par_names = []\n",
    "for i in range(2, 5):\n",
    "    for j in range(4):\n",
    "        x_par_names.append(f'D-{i}T-{j}')\n",
    "        data[f'D-{i}T-{j}'] = data['实际功率'].shift(96*i-j)\n",
    "data = data[96*5:]\n",
    "# x_par_names = ['temperature', 'cloud', 'sun', 'uv', 'water', 'wet']\n",
    "\n",
    "index = data.index\n",
    "columns = data.columns\n",
    "df = mms.fit_transform(data)\n",
    "df = pd.DataFrame(df, index=index, columns=columns)\n",
    "# 预测日为数据集的第365、391天\n",
    "x_train, y_train, x_valid, y_valid = get_train_data_load(df[:365*96], x_par_names, ['实际功率'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LSSVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lssvm import LSSVMRegression\n",
    "\n",
    "# 预测评价指标\n",
    "def get_eval_index(y, y_pre):\n",
    "    # 均方根差RMSE、平均绝对误差MAE、决定系数R2、准确率ACC\n",
    "    rmse = round(metrics.mean_squared_error(y, y_pre)**0.5, 3)\n",
    "    mae = round(metrics.mean_absolute_error(y, y_pre), 3)\n",
    "    r2 = round(metrics.r2_score(y, y_pre), 3)\n",
    "    acc = round(1 - np.mean(np.abs((y - y_pre) / 20)), 3)\n",
    "    return {'rmse':rmse, 'mae':mae, 'r2':r2, 'acc':acc}\n",
    "    \n",
    "# init C = 1, sigma = 1\n",
    "lssvm = LSSVMRegression(c=4.1, sigma=2.8)\n",
    "_ = lssvm.fit(x_train, y_train)\n",
    "y_pre = lssvm.predict(x_valid)\n",
    "\n",
    "_ = mms.fit_transform(load_data[['实际功率']])\n",
    "y_valid_inversed = mms.inverse_transform(y_valid.reshape(-1,1))\n",
    "y_pre_inversed = mms.inverse_transform(np.array(y_pre).reshape(-1,1))\n",
    "_ = plt.plot(y_valid_inversed, 'b')\n",
    "_ = plt.plot(y_pre_inversed, 'r')\n",
    "_ = plt.show()\n",
    "\n",
    "# 预测评价指标\n",
    "eval_index = get_eval_index(y_valid_inversed, y_pre_inversed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 模型对比"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 预测评价指标\n",
    "def get_eval_index(y, y_pre):\n",
    "    # 均方根差RMSE、平均绝对误差MAE、决定系数R2、准确率ACC\n",
    "    rmse = round(metrics.mean_squared_error(y, y_pre)**0.5, 3)\n",
    "    mae = round(metrics.mean_absolute_error(y, y_pre), 3)\n",
    "    r2 = round(metrics.r2_score(y, y_pre), 3)\n",
    "    acc = round(1 - np.mean(np.abs((y - y_pre) / 20)), 3)\n",
    "    return {'rmse':rmse, 'mae':mae, 'r2':r2, 'acc':acc}\n",
    "\n",
    "# RF\n",
    "forest = RandomForestRegressor(\n",
    "  n_estimators=10,\n",
    "  random_state=1,\n",
    "  n_jobs=-1,\n",
    "  max_features='auto',\n",
    "  min_samples_leaf=2, min_samples_split=3)\n",
    "_ = forest.fit(x_train, y_train)\n",
    "forest_y_pre = forest.predict(x_valid)\n",
    "\n",
    "# ANN\n",
    "bpnn = MLPRegressor(\n",
    "  hidden_layer_sizes=(20,),\n",
    "  batch_size=96,\n",
    "  random_state=10, learning_rate_init=0.001)\n",
    "_ = bpnn.fit(x_train, y_train)\n",
    "bpnn_y_pre = bpnn.predict(x_valid)\n",
    "\n",
    "# XGBoost\n",
    "xgboost = XGBRegressor(base_score=0.3, colsample_bylevel=1,\n",
    "    colsample_bytree=0.7, gamma=0, learning_rate=0.05,\n",
    "    max_depth=6, min_child_weight=2, n_estimators=566,\n",
    "    reg_alpha=0.1, reg_lambda=0.05, subsample=0.7)\n",
    "_ = xgboost.fit(x_train, y_train)\n",
    "xgb_resluts = xgboost.predict(x_valid)\n",
    "\n",
    "# MLR\n",
    "linear = LinearRegression()\n",
    "_ = linear.fit(x_train, y_train)\n",
    "linear_y_pre = linear.predict(x_valid)\n",
    "\n",
    "forest_y_pre_inversed = mms.inverse_transform(forest_y_pre.reshape(-1,1))\n",
    "linear_y_pre_inversed = mms.inverse_transform(linear_y_pre.reshape(-1,1))\n",
    "bpnn_y_pre_inversed = mms.inverse_transform(bpnn_y_pre.reshape(-1,1))\n",
    "\n",
    "ts = [0, 1]\n",
    "i = 0\n",
    "t1 = 96*ts[i]\n",
    "t2 = 96*ts[i+1]\n",
    "\n",
    "_ = plt.plot(y_valid_inversed[t1:t2], label='Actual power')\n",
    "_ = plt.plot(forest_y_pre_inversed[t1:t2], label='RF')\n",
    "print(get_eval_index(y_valid_inversed[t1:t2], forest_y_pre_inversed[t1:t2]))\n",
    "_ = plt.plot(linear_y_pre_inversed[t1:t2], label='MLR')\n",
    "print(get_eval_index(y_valid_inversed[t1:t2], linear_y_pre_inversed[t1:t2]))\n",
    "_ = plt.plot(bpnn_y_pre_inversed[t1:t2], label='ANN')\n",
    "print(get_eval_index(y_valid_inversed[t1:t2], bpnn_y_pre_inversed[t1:t2]))\n",
    "_ = plt.plot(y_pre_inversed[t1:t2], label='GAN-LSTM')\n",
    "print(get_eval_index(y_valid_inversed[t1:t2], y_pre_inversed[t1:t2]))\n",
    "_ = plt.legend(fontsize=12)\n",
    "_ = plt.xlabel('time/quarter')\n",
    "_ = plt.ylabel('power/kW')\n",
    "# _ = plt.savefig('/../Image/working-day-result.svg', format='svg')\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = os.listdir(\"./data/多元负荷数据/冷负荷/\")\n",
    "for file in files:\n",
    "    df = pd.read_csv(f\"./data/多元负荷数据/冷负荷/{file}\")\n",
    "    _ = plt.plot(df.iloc[:,1], label=df.columns[1])\n",
    "_ = plt.legend(fontsize=12)\n",
    "_ = plt.xlabel('time/h')\n",
    "_ = plt.ylabel('power/kW')\n",
    "_ = plt.show()\n",
    "\n",
    "files = os.listdir(\"./data/多元负荷数据/热负荷/\")\n",
    "for file in files:\n",
    "    df = pd.read_csv(f\"./data/多元负荷数据/热负荷/{file}\")\n",
    "    _ = plt.plot(df.iloc[:,1], label=df.columns[1])\n",
    "_ = plt.legend(fontsize=12)\n",
    "_ = plt.xlabel('time/h')\n",
    "_ = plt.ylabel('power/kW')\n",
    "_ = plt.show()\n",
    "\n",
    "files = os.listdir(\"./data/多元负荷数据/电负荷/\")\n",
    "for file in files:\n",
    "    df = pd.read_csv(f\"./data/多元负荷数据/电负荷/{file}\")\n",
    "    _ = plt.plot(df.iloc[:,1], label=df.columns[1])\n",
    "_ = plt.legend(fontsize=12)\n",
    "_ = plt.xlabel('time/h')\n",
    "_ = plt.ylabel('power/kW')\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 区间预测"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_interval(arr, dis, conf_level):\n",
    "    '''\n",
    "    arr: array数组\n",
    "    dis: 分布类型, 例如'student', 'normal'\n",
    "    '''\n",
    "    interval = None\n",
    "    mean = arr.mean()\n",
    "    std = arr.std()\n",
    "    if dis == 'student':\n",
    "        interval = stats.t.interval(conf_level, len(arr)-1, mean, std)\n",
    "    elif dis == 'normal':\n",
    "        std = arr.std(ddof=1)\n",
    "        interval = stats.norm.interval(conf_level, mean, std)\n",
    "    return interval\n",
    "    \n",
    "# 预测误差\n",
    "pre_error = y_valid_inversed - y_pre_inversed\n",
    "# 预测区间\n",
    "levels = [0.95, 0.90, 0.85, 0.80, 0.75, 0.70]\n",
    "intervals = []\n",
    "for level in levels:\n",
    "    interval = get_interval(pre_error, 'student', level)\n",
    "    print(f'{level}预测下界为: {interval[0]}, {level}预测上界为:{interval[1]}')\n",
    "    intervals.append(interval)\n",
    "\n",
    "down_data = [[] for i in range(6)]\n",
    "up_data = [[] for i in range(6)]\n",
    "for i in range(96):\n",
    "    down_data[0].append(y_pre_inversed[t1:t2][i][0] + intervals[0][0])\n",
    "    up_data[0].append(y_pre_inversed[t1:t2][i][0] + intervals[0][1])\n",
    "    down_data[1].append(y_pre_inversed[t1:t2][i][0] + intervals[1][0])\n",
    "    up_data[1].append(y_pre_inversed[t1:t2][i][0] + intervals[1][1])\n",
    "    down_data[2].append(y_pre_inversed[t1:t2][i][0] + intervals[2][0])\n",
    "    up_data[2].append(y_pre_inversed[t1:t2][i][0] + intervals[2][1])\n",
    "    down_data[3].append(y_pre_inversed[t1:t2][i][0] + intervals[3][0])\n",
    "    up_data[3].append(y_pre_inversed[t1:t2][i][0] + intervals[3][1])\n",
    "    down_data[4].append(y_pre_inversed[t1:t2][i][0] + intervals[4][0])\n",
    "    up_data[4].append(y_pre_inversed[t1:t2][i][0] + intervals[4][1])\n",
    "    down_data[5].append(y_pre_inversed[t1:t2][i][0] + intervals[5][0])\n",
    "    up_data[5].append(y_pre_inversed[t1:t2][i][0] + intervals[5][1])\n",
    "\n",
    "# _ = plt.figure(figsize=(14, 8))\n",
    "_ = plt.plot(range(96), y_valid_inversed[t1:t2], label='Actual power', color='r')\n",
    "_ = plt.plot(range(96), y_pre_inversed[t1:t2], label='Prediction value', color='black')\n",
    "_ = plt.fill_between(range(96), down_data[0], up_data[0], color='#00ffcc', alpha=0.3, label=f'95%confidence')\n",
    "_ = plt.fill_between(range(96), down_data[1], up_data[1], color='#00cccc', alpha=0.3, label=f'90%confidence')\n",
    "_ = plt.fill_between(range(96), down_data[2], up_data[2], color='#0099cc', alpha=0.3, label=f'85%confidence')\n",
    "_ = plt.fill_between(range(96), down_data[3], up_data[3], color='#0066cc', alpha=0.3, label=f'80%confidence')\n",
    "_ = plt.legend(fontsize=12)\n",
    "_ = plt.xlabel('Time/quarter')\n",
    "_ = plt.ylabel('Power/kW')\n",
    "_ = plt.savefig(work_dir + '/../Image/working-day-interval.svg', format='svg')\n",
    "_ = plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env_b",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "a1d72e55f37f12004b9b2d82a45c015a23cdb2733c870ac1c473b30c47b0950e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
